function fent = nlev_u(nfield,Afield) 
%clear all
np= 4;
NT = 10000;
j = ((np)-1)/2;
m = -j;
jy = zeros(4);
jz = 0;
          for ii=1:np-1
          jz2(ii) = m^2;
          jz(ii,ii) = m;
          jnoz(ii) = sqrt(j*(j+1)-m*(m+1));
          jy(ii,ii+1) = jnoz(ii)*i;
          jy(ii+1,ii) = -jnoz(ii)*i;
          jx(ii,ii+1) = jnoz(ii);
          jx(ii+1,ii) = jnoz(ii);
          m = m+1;
          end  
          jz2(np) = j^2;
          jz(np,np) = j;
          jy=jy/2;
          jx = jx/2;
          rrho0 = zeros(2);
          rrho0(1,1) = 1.;
          sigz = 0.5*([-1 0 ; 0 1]);
          sigx = 0.5*([0 1 ; 1 0]);
          sigy = 0.5*([ 0 i ; -i 0 ]);
          
        
%      parameters
         delta0 = 3.e-3;  % energy difference
         delta1 = 3.e-3; % second bit
        % uz2 = 2*delta0/j;
        %  lam = 1.e1; % field constraint
        %  fieldcon = 0;
%          finding evecs & evals for H0:
          
          E0(1)= -delta0-delta1;
          E0(2)= -delta0+delta1;
          E0(3)= delta0-delta1;
          E0(4) = delta0+delta1;
% field coupling operator
          muop= zeros(4);
          muop(1,1) = 1.;
          muop(4,4) = 1.;
          muop(2,3) = 1.;
          muop(3,2) = 1.;
          %          H0op = trace(delta0*jz*rrho0);
%          L0op = trace(-delta0*jx*rrho0);
%          C0op = trace(jy*rrho0);
        
%         H0 = uz2*diag(jz2)+ delta0*jx;
         
          %for ii=1:np-1
         % H0(ii,ii) = uz2*jz2(ii);
         % H0(ii,ii+1) = delta0*jnoz(ii);
         % H0(ii+1,ii) =  H0(ii,ii+1)';
         % end 
         % H0(np,np) = uz2*jz2(np);
 %        H0 = delta0*jz;
       % call diagu(np,H0,Uop1,E0)
%       [Uop1,E0 ] = eig(H0);
%       E0=diag(E0);
%       jz = Uop1*jz*Uop1';
%       jy = Uop1*jy*Uop1';
       %jx = Uop1*jx*Uop1';
       
       %            jz = matmul(matmul(H0,jz),conjg(transpose(H0)))
%            jy = matmul(matmul(H0,jy),conjg(transpose(H0)))
           
        beta = 1e0/delta0;   % Temperature
        bathcoup = 1e-15;  % e6 for decay to equiv. 
        bpower = 3; % power of omega. depend on the bath model. 3 for the vacuum.
%      light params
        twidth = 1.e-14/2.4e-17 ;
        trun = 2e0*twidth;
        tcenter = trun/2;
        t = trun/(NT-1);
        om0 = 3*2*pi/trun;
        mu = -om0/delta0;
        smu = sqrt(mu^2+1);
        
%      inertial parameters:
        inrtmu0 = -0.5e0;
        inrtalpha0 = pi*1e-8;
        inrtdelta = 0 ; %-inrtalpha0       
          inrtgamma = 0*inrtalpha0; %
 %       load infile
 %       nfield = infile(1);
 %       Afield = infile(2:nfield+1);
%       Afield=0.*Afield;
        U0 = eye(np);
%       U0(2,2) = exp(zi*1d-10) !-1.d0

        nch = 10  ;

%         system-bath coupling operator, currently \sigma_x
         coupop = 0;
         coupop = jy;
%        coupop(1,2) =  1.d0
%        coupop(2,1) =  1.d0

%        time propagation
        time = 0 ;
        l=1 ;
        Uop = U0 ;
           
%       density operator, in the interaction picture:
% eq. relations:
% Pg/Pe = 1+exp(beta*delta-1)
         rho = zeros(4) ;
         rho1 = zeros(4);
         rho2= zeros(4);
         rho3= zeros(4);
         norm= 0 ;
         ent0 = 0 ;
 %         for ii=1:np
        rho(1,1) = 1.;
        rho1(2,2) = 1.;
        rho2(3,3) = 1.;
        rho3(4,4) = 1.;
%        rho(2,2) = exp(-beta*delta0);
 %         rho(ii,ii) = exp(-beta*E0(ii));
 %         norm = norm+rho(ii,ii);
 %         end 
%         rho = rho/norm;
         for ii=1:np
          ent0 = ent0-rho(ii,ii)*log(rho(ii,ii));
         end 
         condswap = zeros(4);
         condswap(1,1) = 1;
         condswap(4,4) = 1;
         condswap(2:3,2:3) = 1/2*( [1+i 1-i ; 1-i 1+i]);
         
          rhof = condswap*rho*condswap';
          rhof1 = condswap*rho1*condswap';
          rhof2 = condswap*rho2*condswap';
          rhof3 = condswap*rho3*condswap';
        emax = 0;
        omega(1:np*np) = 0;
        epsi(1:np) = 0;
        cvls(1:np,1) =1 ;
        while (time<trun)
            fieldeps = 0.;
          for ii=1:nfield
         fieldeps = fieldeps+Afield(ii)*sin(time*ii*pi/trun);
          end 
          fieldeps = fieldeps*exp(-((time-tcenter)/twidth)^2);
          field(l) = fieldeps;
%          if (abs(fieldeps).lt.1.e-80) fieldeps = 0.d0
 %          Hop = fieldom*jz+fieldeps*jx;
 %          Lop = fieldeps*jz-fieldom*jx;
 %          Cop = jy;
 %          theta = delta0*time;
           
            [Uop,emax] = propag2Q(Uop,t,nch,np,fieldeps,emax,E0,muop);
                Uop2 = Uop;
 
                unita(l) = trace(Uop*Uop');
%                expectx(l) = trace(Uop*jx*Uop'*rrho0);
%                expecty(l) = trace(Uop*jy*Uop'*rrho0);
%                expectz(l) = trace(Uop*jz*Uop'*rrho0);
%                expectH(l) = real(trace(Uop'*Hop*Uop*rrho0));
%                expectL(l) = real(trace(Uop'*Lop*Uop*rrho0));
%                expectC(l) = real(trace(Uop'*Cop*Uop*rrho0));
%                Hout(l) = (mu^2*cos(theta*smu)+1)/smu^2*H0op;
%                Hout(l) = Hout(l) + mu*sin(theta*smu)/smu*L0op;
%                Hout(l) = Hout(l) + mu*(1-cos(theta*smu))/smu^2*C0op;
%                Lout(l) = -mu*sin(theta*smu)/smu*H0op;
%                Lout(l) = Lout(l) + cos(theta*smu)*L0op;
%                Lout(l) = Lout(l) + sin(theta*smu)/smu*C0op;
%                Cout(l) = mu*(1-cos(theta*smu))/smu^2*H0op;
%                Cout(l) = Cout(l) - sin(theta*smu)/smu*L0op;
%                Cout(l) = Cout(l) +  (cos(theta*smu)+mu^2)/smu^2*C0op;
%               Uop1 = matmul(Uop1,Uop2)
%               cnorm = trace(2,matmul(Uop1,transpose(conjg(Uop1))))
%            write(17,112) time, real(trace(2,matmul(Uop2,transpose(conjg(Uop2))))),fieldeps
%            write(17,112) real(time,4), real(abs(Uop2(1,1))**2,4),real(abs(Uop1(2,1))**2,4), real(fieldom,4)
%            write(18,112) real(time,4), real(abs(Uop2(2,2))**2,4),real(abs(Uop1(1,2))**2,4), real(fieldeps,4)
%        call diagu(np,Uop2,Uop1,evls)
        [Uop2,evls] = eig(Uop2);
        evls= diag(evls);
        
     
  %           epsi = atan(imag(evls)./(real(evls)));
    
%        arranging the right order of evls and efuncs
           if (l > 1)
            for ii1=1:np
            psi = Uop2(:,ii1);
             spot = 0;
             movlp = 0.;
             for ii=1:np
              psi1 = OldU(:,ii);
             cnorm =  psi'*psi1;
              if (abs(cnorm)>movlp) 
              spot = ii;
              movlp = abs(cnorm);
              end
             end 
             Uop1(:,spot) = Uop2(:,ii1);
             auxevls(spot) = evls(ii1);
            end 
%           Uop  = Uop1
            Uop2 = Uop1;
              
               cvls(:,l) = auxevls;
                    for ii=1:np
     %      epsil(ii,1:l) = phase(cvls(ii,1:l));
     %      epsil(ii,l) = real(i*log(cvls(ii,l)));
             epsil(ii,l) = oldepsi(ii)-real(i*diff(cvls(ii,l-1:l))/cvls(ii,l));
                    end
         epsi = epsil(:,l);
               
            end
            oldepsi = epsi;
            oldomega = omega;
            OldU = Uop2;
            epsiout(:,l) = epsi;
%      write(4000,112) epsi
%         creating the complete set of operators and evls:
            iii=1;
              UT = Uop2';
              for ii1=1:np
                for ii=1:np
                trow = Uop2(1:np,ii);
                tcol = UT(ii1,1:np);
                obasis(:,:,iii) = trow*tcol;
                
%      write(2000+iii,112) real(obasis(:,:,iii))
%      write(3000+iii,112) aimag(obasis(:,:,iii))

               omega(iii) = epsi(ii1)-epsi(ii);
               alpha(iii) = (omega(iii)-oldomega(iii))/t;
               nnumber(iii) = 0.;
               if (abs(alpha(iii))>1d-12)
%              nnumber(iii) = (exp(beta*abs(alpha(iii))-1.d0))**(-1)
               nnumber(iii) = (exp(beta*abs(alpha(iii)))-1)^(-1);
%                 if (alpha(iii).lt.0.) nnumber(iii)=nnumber(iii)+1.d0
                  if (alpha(iii) < 0 ) 
                      nnumber(iii)=nnumber(iii)+1.d0;
                  end
               nnumber(iii) = nnumber(iii)*bathcoup*abs(alpha(iii))^bpower;
               end 
                iii=iii+1;
                end 
                end
     
         
%        checking completness:
%            do i=1,np*np
%               Uop1 = obasis(:,:,i,l)
%               Uop1 = transpose(conjg(Uop1))
%            do ii=1,np*np
%               Uop2 = matmul(obasis(:,:,ii,l),Uop1)
%               comp(i,ii,l) = trace(np,Uop2)
%            end do
%            end do
%              call random_number(rnums)
%            psi = rnums(1:np*np) + zi*rnums(np*np+1:2*np*np)
%                       write(16,*) real(time,4), real(cnorm)
%            call vprod(psi,psi,np*np,cnorm)
%           psi = psi/sqrt(abs(cnorm))
%            call vprod(psi,psi,np*np,cnorm)
%           psi1 = matmul(comp(:,:,l),psi)
%            call vprod(psi1,psi,np*np,cnorm)
             
%          operator expansion coeffcients:
            for ii=1:np*np
            UT = obasis(:,:,ii)';
            coef(ii) = abs(trace(UT*coupop))^2;
            coef(ii) = coef(ii)*nnumber(ii);
            end
            
            gamax = abs(max(coef));
            if abs(gamax<1e-12)
                gamax = emax/1.e4;
            end

          rho = inprop(emax,gamax,np,t,rho,obasis,coef);     
          rho1 = inprop(emax,gamax,np,t,rho1,obasis,coef);
          rho2 = inprop(emax,gamax,np,t,rho2,obasis,coef);
          rho3 = inprop(emax,gamax,np,t,rho3,obasis,coef);    
          
            rhosch = Uop*rho*Uop';
            rhosch1 = Uop*rho1*Uop';
            rhosch2 = Uop*rho2*Uop';
            rhosch3 = Uop*rho3*Uop';
            
            sig1(l,1)=trace(sigx*rhosch1(2:3,2:3));
            sig1(l,2)=trace(sigy*rhosch1(2:3,2:3));
            sig1(l,3)=trace(sigz*rhosch1(2:3,2:3));
            sig2(l,1)=trace(sigx*rhosch2(2:3,2:3));
            sig2(l,2)=trace(sigy*rhosch2(2:3,2:3));
            sig2(l,3)=trace(sigz*rhosch2(2:3,2:3));
          
            
            

%            entropy 
            temprho = rho;
            rhoeig=eig(temprho);
            entropy(l) = -sum(real(rhoeig).*log(real(rhoeig)));
            temprho = rho1;
            rhoeifg=eig(temprho);
            entropy1(l) = -sum(real(rhoeig).*log(real(rhoeig)));
            temprho = rho2;
            rhoeig=eig(temprho);
            entropy2(l) = -sum(real(rhoeig).*log(real(rhoeig)));
            temprho = rho3;
            rhoeig=eig(temprho);
            entropy3(l) = -sum(real(rhoeig).*log(real(rhoeig)));
            
            
            for iii=1:np
             pop(iii,l) = real(rhosch(iii,iii));
             pop1(iii,l) = real(rhosch1(iii,iii));
             pop2(iii,l) = real(rhosch2(iii,iii));
             pop3(iii,l) = real(rhosch3(iii,iii));
            end 
            fid(l) = real(trace(rhosch'*rhof));
            fid(l) = fid(l) + real(trace(rhosch1'*rhof1));
            fid(l) = fid(l) + real(trace(rhosch2'*rhof2));
            fid(l) = fid(l) + real(trace(rhosch3'*rhof3));
        
    %        fieldcon = fieldcon+lam*fieldeps^2;

        time=time+t ;
        l=l+1;
              end 
     4-fid(l-1);
     fent = fid(l-1);
%     1-fent ;    
 %          write(666,112) entropy+fieldcon
